Phase transitions for rock-scissors-paper game on different networks 
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Monte Carlo simulations and dynamical mean-field approximations are performed to study the phase tran- 
sitions in rock-scissors-paper game on different host networks. These graphs are originated from lattices by 
introducing quenched and annealed randomness simultaneously. In the resulting phase diagrams three different 
stationary states are identified for all structures. The comparison of results on different networks suggests that 
the value of clustering coefficient plays an irrelevant role in the emergence of a global oscillating phase. The 
critical behavior of phase transitions seems to be universal and can be described by the same exponents. 
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Complex networks are recently studied extensively as these 
objects pervade all nature Ql y[ yj. Initial efforts were fo- 
cused on the characterization and evolution of these graphs. 
Very recently the investigations are extended to critical behav- 
ior in dynamical models defined on various random networks 
|4j,|5|]. The effects of the network topology on critical tran- 
sitions can well be investigated on the small-world networks 
introduced by Watts and Strogatz [6] because this structure 
provides a transition from a regular lattice to some random 
networks. It is shown that the modifications in the host topol- 
ogy either change the class of universality 013] or the long- 
range connections prevent complete ordering [9]. 

Besides it, there are examples where the increase of ran- 
domness of the host network induces a phase transition to- 
wards a state that cannot be observed on lattices. To be 
specific, the emergence of an oscillating state has been ob- 
served by Kuperman and Abramson llOll at a finite rate of 
randomness in an epidemiological model. The appearance of 
global oscillation (synchronization) is also observed in some 
cyclically dominated three-state systems, such as the volun- 
tary prisoner's dilemma full and rock-scissors-papers (RSP) 
games 1 12]. For the latter models the amplitude of oscillation 
depends on the randomness and a second transition may occur 
when the limit cycle approaches the absorbing states. The ap- 
plication of annealed (temporal) randomness results in quali- 
tatively similar behaviors. In this Brief Report, concentrating 
on RSP game, we extend our previous work by combining 
both types of randomness to explore the global phase diagram 
and to clarify a universal feature. 

For the spatial rock-scissors-paper game the individuals lo- 
cated on site i of a lattice belong to one of the three species 
(sj = 1,2,3) which dominate cyclically each other. This 
means that the time evolution of this system is governed by 
the iteration of cyclic invasion processes between two ran- 
domly chosen neighboring sites. Namely, the pairs (1, 2) or 
(2, 1) become (1, 1), (2, 3) or (3, 2) become (2, 2), and finally 
(3, 1) or (1, 3) evolve to (3, 3) with the same rate defined to 
be 1 . Starting from a random initial state on a square lattice 
this system evolves into a self-organizing pattern in which 
the three species are present with the same average concen- 
tration (1/ 3). This state maintained by the cyclic invasions 
I13| 14 lll5ll provides protection against some types of invaders 
flrj[l7ll and exhibits an unusual response to the external per- 
turbation imm. 



In our previous paper 1 12] the RSP game is studied on 
such regular small-world structures where randomly chosen 
(long-range) links are substituted for portion Q of the nearest- 
neighbor links between sites located on a square lattice. A 
similar algorithm was proposed by Watts and Strogatz jg] as a 
continuous transition between a lattice and a random network. 
For the sake of simplicity, our rewiring process 11211 conserves 
the regularity, i.e. the coordination number of each site (called 
degree) remains unchanged (z = 4). Evidently, for Q = this 
structure is equivalent to a square lattice meanwhile the limit 
Q — > 1 results in a regular random graph. It is important to 
note that the restriction of fixed degree does not change the 
small-world feature of the network: the characteristic (aver- 
age) path length between two sites scales as the logarithm of 
the network size. Consequently the main conclusion is ex- 
pected to be also valid for the standard Watts-Strogatz net- 
works. At the same time, this simplification has allowed us 
to compare the results from these quenched structures with 
those obtained from "annealed" structures that can also be in- 
vestigated analytically. For the annealed structures the (long- 
range) links are substituted temporarily with a probability P 
for the standard (nearest-neighbor) links on the square lattice. 
Naturally, in the limit P — ► 1 the dynamics of the system 
satisfies the mean-field condition. The effect of annealed ran- 
domness on rumor propagation 1 20, 21] and for the two-state 
voter model 1 22] has already been studied. 

For the RSP model on the quenched structures the emer- 
gence of global oscillation is observed if Q exceeds a thresh- 
old value and the amplitude of oscillation tends to a fixed 
value in the limit Q — ► 1. Such a transition occurs also for 
the annealed structures when increasing the value of P. In the 
latter case, however, the amplitude of global oscillation in- 
creases with time and finally the evolution ends in one of the 
homogeneous (absorbing) states if P exceeds a second thresh- 
old value. The different consequences of the two types of ran- 
domness have inspired us to combine these randomness when 
starting from different lattices. 

Henceforth, both types of randomness are applied to derive 
a global phase diagram on the P — Q plane. To explore the 
robust (universal) behaviors we also analyze what happens for 
other regular structures with degree of z = 3, 4, and 6. These 
lattices are Kagome (where z = 4 as for the square lattice), 
honeycomb (z — 3), triangle (z = 6), cubic (z = 6), and 
ladder-shape (z — 3) lattices. For the ladder-shape structure 
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FIG. 1 : Phase diagram for the RSP game on network derivated from 
honeycomb lattice. The system exhibits self-organizing (S), oscillat- 
ing (O), and absorbing (A) phases. Symbols represent MC results. 
The solid lines are guides to the eye. 



two parallel chains are connected by interchain bonds. 

Our analysis is based on systematic Monte Carlo (MC) sim- 
ulations. First we create a quenched random regular structure 
starting from one of the above mentioned lattices. That is, a 
portion Q of the nearest-neighbor links are replaced by ran- 
domly chosen links in a way that conserves the regularity (for 
details see O). The MC simulations are started from a ran- 
dom (uncorrelated) initial state where the three species take 
their place with the same probability (1/3). Keeping this struc- 
ture fixed the time evolution is governed by invasions between 
neighbors with a probability 1 — P or along a (random) long 
range link chosen with a probability P. In the simulation the 
number of lattice points is varied from 10 3 to 10 7 . The large 
sizes are used in the close vicinity of transition points to re- 
duce the undesired effect of fluctuations. 

The main feature of the steady-state phase diagram, which 
is generally valid for all structures, can be summarized as fol- 
lows. For small values of Q and P, the stationary-state is char- 
acterized by a self-organizing strategy distribution denoted by 
S. In this self-organizing pattern the strategies alternate cycli- 
cally at each site, these local oscillations are not synchronized 
by the short range interactions, and the average concentra- 
tions are the same (1/3). For the opposite limit - when both 
structural parameters, Q and P, are close to 1 - the system 
evolution is characterized by growing spiral trajectories I12ll 
and finally the evolution ends in one of the three absorbing 
(homogeneous) states (A) containing only one strategy. Evi- 
dently, this absorbing phase is three-fold degenerated due to 
the cyclic symmetry. These two phases (S and A) are sep- 
arated by the region of global oscillation (O) on the P — Q 
plane. In this oscillating phase the behavior is characterized 
by a limit cycle which is quantified by an order parameter $ 
defined as the ratio of the area of the limit cycle and area of 
the full triangle in the ternary phase diagram 1 12]. This order 
parameter is (1) in the state S (A) and varies from to 1 
for the occurrence of global oscillation (phase O). Figure ^ 
shows the phase boundaries obtained by MC simulations for 
2 = 3 on the P — Q plane. 



FIG. 2: The order parameter as a function of Q for quenched net- 
works (P = 0) derived from cubic (box), square (bullet), and ladder- 
shape lattice (triangle). The corresponding critical value of Q where 
the oscillatory phase emerges is Qi (P = 0) = 0.0103(1), 0.067(1), 
and 0.210(2) respectively. The curves from left to right in the inset 
show the $((3, P = 0) functions for small-world networks derived 
from cubic, triangle, square, Kagome, honeycomb, and ladder-shape 
lattices. 



A striking quantitative difference is found in the behavior of 
these systems when varying the value of z. Namely, global os- 
cillation (phase O) can be observed on random regular graphs 
(Q = 1, P = 0) for z = 3 and z = 4. On the contrary, the 
system evolution terminates in one of the absorbing states (A) 
on random regular graphs for z — 6 and similar behavior is 
expected for z > 6. In order to illustrate the relevance of z 
Fig. |21 shows the Q-dependence of <I> for quenched structures 
(P = 0) created by the rewiring technique from the above 
mentioned lattices. 

Evidently, in the limit Q — > 1 the quenched structures be- 
come independent of the original lattice. Consequently, start- 
ing from either the square or the Kagome lattices the order pa- 
rameter $((3, P = Q) tends monotonously to the same value 
($i(z = 4) = 0.980(5)) in the limit Q -> 1. For z = 3 (trian- 
gular and ladder-shape lattices) we have obtained a lower limit 
value, namely <&i(z = 3) = 0.750(5). For z = 6, however, 
the system ends in the phase A if Q exceeds a threshold value 
Q2- According to our simulations on the quenched structures 
(P = 0) Q2 = 0.378(8) if the random regular structure is 
created from a cubic lattice and Q2 = 0.405(10) for the trian- 
gular lattice. 

The topological structure of the original lattice affects 
the value of Qi because the global oscillation occurs for 
small Q for all the investigated lattices. The inset in Fig. |2 
demonstrates clearly that the lowest Q\ appears for the three- 
dimensional cubic lattice, while the highest value is found for 
the one-dimensional ladder-shape structures. It is underlined 
that the values of Q\ are very close to each other for all the 
small-world structures created from the two-dimensional lat- 
tices. 

The present RSP system undergoes two subsequent phase 
transitions when increasing the randomness (Q and/or P) of 
the backgrounds. The first transition (from S to O) is a Hopf- 
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bifurcation that is well studied by mean-field type approaches 
ll23ll . Notice that the order parameter $ vanishes linearly for 
all the structures (see Fig.|2j in agreement with the theory. 

Our results indicate clearly that the global oscillation 
emerges just above a threshold value of quenched randomness 
although the small-world feature characterizes this type of 
networks at any small rate of disorder. The same phenomenon 
was observed by Kuperman and Abramson 1 10] when con- 
sidering the transition to an oscillating phase in a three- 
state (susceptible-infected-refractory) epidemiological model 
where the initial lattice was the traditional one-dimensional 
ring. They have conjectured that the emergence of global os- 
cillation is related to the variation in the clusterization. In 
the present model, however, our results suggest that the clus- 
tering coefficient (C) cannot play significant role because the 
corresponding values of C are very different for the studied 
structures. We cite as an example the case where C = for all 
initial lattices except the triangle and Kagome where C = 2/5 
and 1/3 respectively. Furthermore, the value of C is ranging 
from 0.0004 to 0.3 for different topologies when the system 
enters into the global oscillating phase. The inset of Fig. |2] 
shows the order parameter as a function of Q for all networks 
studied here. These results suggest that Q\ depends basically 
on the dimension of the initial lattice. Namely, Q\ = 0.068(6) 
for all networks derived from a two-dimensional lattice inde- 
pendently of the value of coordination number. On the con- 
trary, Qi is much smaller for the graph derivated from the cu- 
bic lattice and substantially larger for the network originated 
from the ladder-shape (practically one-dimensional) lattice. 

The second phase transition (from the state O to ^4) can be 
studied more efficiently (or with a higher accuracy) if vary- 
ing P rather than Q because the limit cycle is affected by the 
quenched structural randomness even for large N. This dis- 
crepancy can be avoided by averaging over many runs on dif- 
ferent structures whose creation is very time-consuming. For 
temporary randomness (P > 0), however, the technical diffi- 
culties can be overcome more easily and the numerical analy- 
sis is executable with an adequate accuracy. This is the reason 
why henceforth we concentrate on the second phase transition 
occurring with the increase of P. 

First we study what happens on random regular graphs 
(limit Q — > 1) for z — 3 when increasing P. This choice was 
motivated by the simplicity of this tree-like structure on which 
the dynamical cluster technique can be applied for sufficiently 
large clusters. In Figure[3]the MC results show that the order 
parameter $ tends very slowly to 1. Despite the large sizes 
in simulations we couldn't study the very close vicinity of the 
transition point (P2) because of the fluctuations yielding an 
occasional transition from the noisy limit cycle to one of the 
homogeneous states (^4). Our MC data can be well approxi- 
mated by a power law behavior, namely 1 — $ oc (P2 — P) 7 
with an exponent 7 = 3.3(4). 

For periodic structures the dynamical cluster technique is 
proved to be a very efficient method to describe different phe- 
nomena in several non-equilibrium models I IT2L l24l l25l l2al . 
In the limit TV — > 00 the random regular graph become locally 
tree-like and it can considered as a Bethe lattice on which this 
technique works well too 1 27] . Using this technique one can 
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FIG. 3: The order parameter as a function of P for the three-degree 
random regular graph. Boxes indicate the results of MC simulation. 
The estimated critical value of the transition to absorbing state is 
P2 = 0.079(2). The dashed (solid) line represents the prediction of 
dynamical cluster technique at the levels of four- and six-site approx- 
imation. Inset shows the shape of four- and six-site clusters used in 
these approximations. 



determine the probability of each configuration occurring on 
a fe-site cluster by solving a suitable set of master equations 
(details are given in 

dill). The cluster size k is a crucial 
parameter. For one-site clusters this technique is equivalent 
to the mean-field approximation predicting concentric orbits 
independent of P 12311 - The pair approximation (k = 2) gives 
spiral trajectories reaching the edge of triangles (or the ab- 
sorbing states) for arbitrary P. Choosing four- and six-site 
clusters (see the inset in Fig. [5} this method was capable to re- 
produce the appearance of limit cycles below a threshold value 
of P. The quantitative predictions of this method are com- 
pared with the MC results in Fig. [3] Obviously, the increase 
of cluster size improves the estimation. Notice that according 
to the six-site approximation the order parameter $ also tends 
to 1 very slowly (see, Fig.0}. The extrapolated critical values 
are P 2 (4s) = 0.019(2) and P 2 (6s) = 0.067(4). 

The above value of the exponent 7 agrees (within statisti- 
cal error) with those we have observed on the square lattice in 
our previous work concentrated only on the effect of annealed 
randomness. This coincidence inspired us to study the robust- 
ness of this transition. For this purpose the MC analyses of the 
second transition were carried out on some two-dimensional 
lattices (square, triangle, and honeycomb for Q = 0). The 
results, summarized in Fig.|4] seem to confirm that the tran- 
sition from the global oscillation (phase O) to the absorbing 
states (A) is universal. 

To summarize, we have studied the effect of host lattice 
randomness on the stationary state for a simple rock-scissors- 
paper system. The quenched and annealed randomness of 
the regular background is characterized by two parameters Q 
and P varying from (corresponding to a lattice with a di- 
mension of 1, or 2, or 3) to 1 (random regular graph and/or 
mean-field condition). This system displays two subsequent 
transitions if the randomness is increased. A self-organizing 
pattern can be observed if these randomness parameters re- 
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FIG. 4: Comparison of the continuous transitions to absorbing phase 
for different networks. The MC data are obtained on honeycomb 
(triangle), square (cross), triangle (bullet) lattices, and on random 
regular graph with a degree of z — 3 (box). Filled triangle represents 
the result of dynamical cluster approximation at six-site level for the 
same graph. The solid line represents the slope of f3 = 3.3. 



"amplitude") increases with the randomness up to its satura- 
tion value. Thus, above a second threshold value [more pre- 
cisely for P > PziQ)] the system sooner or later terminates 
in one of the homogeneous absorbing states. Our simulations 
indicate that this second transition has also a universal fea- 
ture in the slow tendency towards the saturation value. It also 
turned out that the global oscillation is stable ($ < 1) on the 
quenched random regular structures if the number of neigh- 
bors is not larger than four (z < 4). It would be interesting 
to see how the synchronization (as well as the above men- 
tioned two transitions) emerges on other random networks. 
The present rock-scissors-paper model involves two crucial 
features. On the one hand, the dynamical rule is cyclically 
symmetric; on the other hand, the invasion fronts become very 
irregular even on the two-dimensional lattices because the in- 
vasion between two neighboring sites is not affected by their 
neighborhood. Further investigations are required to clarify 
what happens if the dynamical rules are not cyclically sym- 
metric and/or the moving interfaces are smoothed by local in- 
teractions. 



main within a region of P — Q plane. When crossing the 
boundary of this region a global oscillation (limit cycle) oc- 
curs via a Hopf-bifurcation. The transition point (for P = 0) 
depends strongly on the dimension d of the original lattice 
meanwhile it is hardly affected by the clustering coefficient. 
For the global oscillation the area of limit cycle (as well as the 
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